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Abstract 

Presented here is a study of well-posedness and asymptotic stability of a “degenerately 
damped” PDE modeling a vibrating elastic string. The coefficient of the damping may vanish 
at small amplitudes thus weakening the effect of the dissipation. It is shown that the resulting 
dynamical system has strictly monotonically decreasing energy and uniformly decaying lower- 
order norms, however, is not uniformly stable on the associated finite-energy space. These 
theoretical findings were motivated by numerical simulations of this model using a finite element 
scheme and successive approximations. A description of the numerical approach and sample 
plots of energy decay are supplied. In addition, for certain initial data the solution can be 
determined in closed form up to a dissipative nonlinear ordinary differential equation. Such 
solutions can be used to assess the accuracy of the numerical examples. 


1 Introduction 

Advances in nonlinear functional analysis and the rich theory of linear distributed parameter sys¬ 
tems have led to a growing of body of work on nonlinear infinite-dimensional models. For instance, 
in a 2nd-order evolution framework (especially, wave, elastodynamics, or thin plates with no rota¬ 
tional inertia terms) for an appropriate elliptic operator A a linear equation with viscous damping 
for an unknown u = u{t, x) may be expressed as 

ii + A{x)u + /3{x)u = F{t, x) 

with /3 > 0. We will focus on the evolution on a bounded domain and under suitable homogeneous 
boundary conditions. A nonlinear refinement on the dissipative term may take the form of a 
feedback law g{u). Stability properties of such models have been extensively analyzed. In an 
infinite-dimensional sitting such a nonlinear feedback may change the topology of the problem and 
uniform stability becomes reliant on the regularity of solutions (for example, see [HEllS]). 
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A more general scenario would account for coefficients that depend on the solution itself: 

u + A{x)u + j3{x, u, u)u = F{t, x). 

Assuming the well-posedness of an associated initial-boundary value problem can be resolved, if the 
term j3{x,u^u) is not guaranteed to be strictly positive on a fixed appropriately configured set, then 
analysis of stability becomes much more involved since the region where the dissipation is active 
now evolves with the solution and may not always comply with the requirements of the geometric 
optics. The case when this coefficient vanishes at zero displacement, namely, /3(x, 0, ii) = 0, will be 
referred to as degenerate damping. 

Such a degeneracy naturally arises when investigating energy decay of higher-order norms. For 
example, the natural energy space for a semilinear wave problem 


u — Au + g{u) = 0 


is u G IT^’^(n) and ii G Lp‘{Q). With more regular initial data one can consider behavior of 
higher-order energy norms, namely for {u,ut) G x VF^’^(n). One approach would be to 

differentiate the PDF in time which via the substitution v = ii leads to a degenerately damped 
problem 

V — Av + g {v)v = 0 

A particular example can be observed in the relation between Maxwell’s system and the (vectorial) 
wave equation. For a given medium denote the electric permittivity by e, magnetic permeability 
by A and conductivity by a. Then Maxwell’s system reads 

E — cml{XH) -|- aE = 0 
H -|- curl(eF) = 0 , 

with div{XH) = 0. On a bounded domain, subject to the eleetrie wall boundary conditions, and for 
scalar-valued A, e, a with positive lower-bounds, the term aE exponentially stabilizes this system 
[3]. In a more accurate nonlinear conduction model the coefficient a may depend on the intensity 
of the electric field E. If we consider, for example, a = for p > 1, then differentiating the 

hrst equation in time and combining with the equation for H gives 

E + curl(Acurl eE) + ap\E\P-^E • EE + a\E\PE = 0 . 

For example, taking, p = 1 gives 

E -|- curl(A curl eE) -|- a{E ■ E)E + a\E\E = i) 

where E is the normalized vector E. The term a\E\E has features of the viscous dissipation in this 
second-order equation, but nonlinear conductivity augments it with a degenerate coefficient a\E\. 

The study of stability for the above models is much more delicate than in the situations where 
the damping, even nonlinear, depends on the time-derivative only. Weighted energy methods—from 
basic energy laws to Carleman estimates (e.g. [a m El El El ini [n]) —have been successfully used 
to derive stabilization and observability inequalities for distributed parameter systems. However, 
these methods typically rely on the properties of the coefficients to ensure that suitable geometric 
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optics conditions are satisfied and the control effect suitably “propagated” [l2] across the physical 
domain. One can sometimes dispense with geometric optics requirements for smooth enough initial 
data da H [El [Ml [IT], yet even then the support of the control/damping term must contain a 
subset that is time-invariant (and with any time-dependent coefficients being non-vanishing, e.g., as 
in m)- In turn, the analysis of control-effect propagation when the coefficients themselves depend 
on the solution and possibly go to zero wherever and whenever the solution does would require new 
techniques. 

1.1 The model 

The following semilinear model, if recast in a higher-dimensional setting becomes highly non-trivial 
even when just regarding local wellposedness. In a one-dimensional framework the nonlinearity 
is more tractable, but the rigorous stability analysis has long been open. We focus on an elastic 
string with a degenerate damping, namely a dissipative term whose coefficient depends on and may 
vanish with the amplitudes; 


u — Uxx + f{u)u = 0, for X £ Q := (0,1), t > 0 (1) 

fixed at the end-points 

tt(t, 0) = u(t, I) = 0 for t>0 (2) 

and with a prescribed initial configuration at t = 0: 

u{0, x) = uo{x), u{0,x) = ui{x) for a.e. x G (0,1). (3) 

The initial data uo,ui live in the natural function spaces revisited below. Function / is assumed 
to be continuous non-negative, hence the term f{u)u a priori should provide some form of energy 
dissipation in the model. 

The scenario of interest is when f{s) —)• 0 as s —)• 0, essentially causing the dissipative effect to 
deteriorate at small amplitudes. We will focus on the polynomial case 

/(s) := a > 0, m G N. (4) 

satisfying the locally Lipschitz estimate 

2m—1 

f{s) — f{r) = M{s,r)-{s — r) with M{s,r) = a ^ . (5) 

j=0 


1.2 Known resnlts and new challenges 

Existence and uniqueness of weak finite-energy solutions to was proven in m by means 

of Galerkin approximations. The advantage of a ID framework is that the displacement function 
is absolutely continuous, hence topologically f{u)u is still in L^(0,1), as in the case of the corre¬ 
sponding linear model. However, in higher-dimensional analogues this embedding property is lost 
and proving existence becomes a markedly more complex task. First, fractional damping exponents 
were considered in order to ensure that the damping term is bounded with respect to the finite 
energy topology |19j . Arbitrary damping exponents were subsequently examined in [201121]. Due to 
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the loss of regularity solutions had to be characterized via a variational inequality and established 
by a rather technical application of Kakutani’s fixed point theorem. 

On the other hand, stability analysis even in one dimension poses a challenge that has been 
open for a number of years. Despite the gain in regularity, attributed to Sobolev embeddings, the 
key difficulty now is that energy estimates require some sort of information on the region where 
the damping is supported. In Q both the magnitude and the support of the damping coefficient 
evolve with the geometry of the state, rendering all standard techniques inapplicable. 

It is plausible to assume that some sort of a logarithmic uniform decay rate can be verified, 
possibly by combining ODE techniques (e.g. [221 I23| l with pointwise Carleman-type estimates. 
Another, though a rather weak, tentative indication of this outcome would be the uniform stability 
of the corresponding finite-dimensional analogue (see the Appendix). Yet the situation in inhnite 
dimensions turns out rather different. 


1.3 Contribution of this work 

The goal of this article is to examine analytically and numerically stability properties of the dy¬ 
namical system associated to 0-®): 

• Establish global persistence of regularity in solutions with smooth initial data. Besides the¬ 
oretical interest such a result is useful to justify the convergence estimates for numerical 
approximations. 

• Prove that a polynomial degeneracy in the damping of the form Q yields a system that is 

not uniformly stable. 

• Present a numerical scheme that indicates the loss in decay rates. Such observations had 
been performed first and, in fact, served as a motivation for the theoretical resnlts presented 
here. 


1.4 Outline 


The notation employed throughout the paper is summarized in Section 2.1 The two main results 


on well-posedness and stability are stated in Section 3.2 


Several auxiliary technical definitions used in the proofs can be found in Section 3.1 Local 


and global wellposedness are verified respectively in Sections |4.4| and 4.5 
regularity lemmas proved earlier in Section [4.3[ 


They draw upon two 


The proof of the lack of uniform stability is the snbject of Section Numerical results are the 
subject of Section 

The Appendix contains results pertaining to the ODE analog of the considered problem, namely, 
a damped harmonic oscillator with the damping coefficient dependent on the displacement. 


2 Preliminaries 

2.1 Notation 

This section serves as a quick reference for the basic notation used thorough the paper with some 
of the symbols revisited and discussed in more detail later in the text. 
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Henceforth || • ||x will denote the norm on a normed space X. For the space L^(n) we will use 

l^lo — II^IIl^(o)) 

with the corresponding inner product denoted by (•, •)(,. We will also frequently involve the Sobolev 
space 

associated to an equivalent inner product and norm 

iu,v)^ = {ux,Vx)o |tt|i := 

The bilinear form (•, •) will indicate the pairing of Hq{Q,) and its continuous dual 
We will also frequently use spaces of the form 

C”([0,r];X) or 1^(0,T;X), 

which will be abbreviated respectively as 

C^X and X. 

Looking ahead, for the one-dimensional Dirichlet Laplacian operator A (discussed below) let us 
introduce the space 

5^ := |z I z G C^j,{3!{A *-2 ^)) for j = 0,1,..., n-|-1)| (6) 

equipped with the natural graph norm. For example, = C([0, T]; n C^([0, T]; L^(n)) 

indicates the standard regularity in space and time for a weak solution to a linear wave equation. 
In turn 5^ = CtSi{A) n n C^L‘^{Q) does the same for a strong solution to such an 

equation. 

We will also be using spaces 

j^n X &{A^). 

Thus is the natural finite energy state space for a linear wave problem and denotes the 
domain of the corresponding evolution generator. 

Relation a < b will occasionally be used to indicate that a < Cb for a constant C which only 
depends on the main parameters of the system, e.g., the size of the domain hi or the exponent p of 
the coefficient in Q. 

2.2 Laplace operator 

For convenience let us summarize some of the fundamentals. Consider the operator 

A =-da,x ■■ ^(A) c L‘^{n) ^ ^{A) = w'^'\n)nH^{n) (7) 

which is positive, self-adjoint with compact resolvent, and has eigenvalues 

A„ = n^vr^, n G N 
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with the corresponding eigenfunctions 

En{x) = \/2sin(n7rx), n G N. 

The EnS form an orthonormal basis for LP‘{VL). For r G M we can define fractional powers of A: 





OO 

n=l 


with 

{ OO OO 

Y^CnEn : 

n=l n=l 

The eigenfunctions {En} form an orthogonal basis for every ^(A''), r G M. Some of the fractional 
powers can be identified with Sobolev spaces, e.g. 



Since in our situation the model is one-dimensional, then trivially no issues in these identifica¬ 
tions arise in regard to the regularity of the domain. Operator A also corresponds to the Riesz 
isomorphism Hq(Q) —)• and for u,v G Hq{Q) we have 

{Au,v) = (u,n)i,. 


2.3 State space and linear group generator 

The natural finite-energy state space associated with the evolution driven by ©-(HI) is 

If we set y = (u, ii) we can recast this problem as an evolution equation 

y = Ay 

for the skew-adjoint operator A : ^(A) C —)• 

A-[° ^1 

with 

^(A) = &{A) X ^ 

We will also consider smoother solutions for which we define 

■.= = &{A^) X &{A^I‘^) for nGN 

with the associated graph norm given via 

11(^0, Ul)||^n =||^'-0||^(^(n+l)/2) + ||^l||^(^n/2) 

OO 

k=l 


( 8 ) 

(9) 
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Here J-k denotes the k-ih Fourier coefficient with respect to the Hilbert basis of L^(H) given by the 
eigenfunctions of A. Note also that in one-dimension the following continuous injection holds 

^ X C”(n) X (10) 

for n G N U {0} if we adopt the notation convention := Lp‘{Q). 

3 Results 

We start with a formalized notion of a solution to the class of PDF systems of the form Q-Q. 
The following slightly more abstract formulation will help streamline the subsequent discussion. 

3.1 Auxiliary definitions 

Definition 3.1 (Wave problem). Let [[T, (/>o, </>i| be the shorthand for the initial-boundary value 
problem: 

4>xx = F, for X € Q, t > 0 , 

with the indicated derivatives taken in the sense of distributions, and subject to boundary conditions 

(j){t, 0) = 0 = (j){t, 1) for t > 0 

and initial data 

x) = 4>o{x), (^(0, x) = (/>i(x) for a.e. x G Ll. 

Definition 3.2 (Weak solution of linear problem). Suppose (uo,ui) G and for some T > 0, 
F G Then we say a function 

u G C''fL^(H) n CtHq{LI) [equivalently u G or [u, ii) G (11) 

is a weak solution to [[F, uo,ui| on interval [0,T] if 

(i) u[i),x) = uo{x) and m( 0,x) = ni(x) for a.e. x G 

(ii) For any (f G LP‘[Vl) the sealar map t i—)■ [u{t),(p)^ is absolutely continuous (hence a.e. differ¬ 
entiable) on [0, T]. 

(Hi) For any cj) G Hq{LI) 

^[u[t),(l))o + {u{t),(j))^ = [F{t),(lf)^ for a.e. t G [0,T). (12) 

Definition 3.3 (“Regular” functions). A function u on [0,T] x H will be described as regular of 
order n G N U {0} if it is continuously differentiable in time with the following regularity 

{u,u)gCt^^. (13) 

In classieal terminology, weak solutions correspond to order 0 and strong solutions to order 1. 


7 


Suppose u = u(i, x) has the weak regularity 0) (regular of order 0). Then according to the 
(ID) Sobolev embedding ^ C'"*“^(D) for m G N, the function F{t,x) = f{u{t,x))u{t,x) 

is well-defined as an element of L|,L^(D). In fact we will generalize this statement for the purposes 
of subsequently analyzing more regular solutions. 

Proposition 3.1. Let f{s) = for m G N. If z € SJf, then 


(14) 


(15) 


A^I’^\f{z)z\^Cij,Ql{A 2 ^)^ for s,jGNU{0}, s-|-j<n. 

In addition, 

\\dtifiz)z)\\crpL2{n)<'P Wdi^WcrHAn) . + \\9tz\\crpL‘2{n)), 3 

L J /c——1 

where V is a polynomial in j — 1 variables. 

Example 3.1. Due to a variety of spaces and indices involved in the statement of Proposition \3A\ 
it is helpful to look at a basic example. Take a = m = 1, so f{z)z = z^z and consider the regularity 
order n = 4. Then the condition z € Stf reads 


z€f] C^^i^iA^)) 

j=0 


(16) 


In particular, z G CT^iA^/"^), which corresponds to 5 square-integrable derivatives, first three of 
which satisfy zero boundary conditions. We have, for example. 


df[f{z)z] = 30iz^ -|- 20z‘^'z + 20z z "z -1- IQzzdfz + z^ 


dfz 


Thus, for instance, \dff{z)z\o can be estimated using a polynomial of L°°{Q) bounds on the functions 
z, z, ..., dfz, and one term involving the Lp‘{Q) norm of the fifth derivative of z in time or, 
equivalently, \dfz\o. This is precisely the conclusion of (15). 

Likewise, if we consider, say, the 2-nd derivative in space and 2-nd in time to s = 2, j = 2 in 
(14) we get: 


df[fiz)z] = 2z^ 6zzz + z^^'z 


d‘f,df[f{z)z] =12zfz -|- 12zxZZx -|- 12zxZZx + ‘Iz'i'z -\- Az'z'xZx -|- 12zZxZx + 6zZxxZ 4- 


1.^2 ■ 


i2 ■ 


-|- 2z'zZxx + (iZZZxx + QzZZxx + Z,^ "z 


(17) 

(18) 


We have that df[f{z)z] has zero trace as follows from and the zero boundary condition on z. 
Moreover, the highest-order term 'z’xx in dxdf[f{z)z] can be bounded in Lp‘{Lt) since by (16) we have 
z G S!{A). The rest of the terms are in fact in L°°(Ll). This confirms that Adf[f{z)z] G CtL‘^{Q.) 
in agreement with (14). 

Note that if we had, say, n = 6 then in the same context we would need to prove that Adf[f{z)z] G 
Ct^{A) instead of just CtL‘^{LI). First of all, n = 6 would give z G CtI^{A’^/'^) implying that 
z, Zx,..., dxZ have zero traces. Hence so do their time-derivatives and then it is immediately fol¬ 
lows that Ad‘f[f{z)z] satisfies zero trace condition. And taking two more derivatives in space in 
(18) yields the highest-order term dfdfz which is in L^(D) from z £ that implies 'z G 
Thus, dfdfz G L^(D) and whereas other summands in dxAdf[f{z)z] belong to L°°{Ll). We conclude 
that ifn = 6, then Adf[f{z)z] (is in L?‘{Ll)) has zero trace and integrable second derivative, again 
in accordance with Q- 













Proof. Let 


Fiz) ■.= f{z)z 

and first let’s show that d^F belongs to For j < n, d^F is a polynomial in the 


following variables 


diF = V[f^^\z), dfz. d!fz]k=ox...„ 


+1 


that is affine with respect to the highest-order derivative d^z = d^^^z, which is at least in CTL'^{Pt) 
by the assertion that z £ SJf (recall Q and plug j = n + 1). We can just bound coefficients of 


z using the fact that -ffo(ri) embeds continuously into C{Q). Thus the bound (15) follows. 


^ _ ... 

Next, S>{A ^2 ^) embeds into for j < n, so all of the terms in dx ^d^F are CTC{Pt) 

except for the highest order term dx~^9lz G CtL‘^{Q). Since dx~^dlF is affine with respect to that 
term with continuous coefficients, then 


It follows that 


d^-^dlF G CtL^{Q) 


d{F G CTW'^-^’^in) for j < 


n. 


To strengthen this regularity to Ct^{A 2 ) we must verify the boundary conditions. Since 
^(A) coincides with the functions that have zero trace, then it is sufficient to show that 


d^djF = 0 on dPt for 0 < A: + j < n — 2 {k,j>0) 


(19) 


That is, we can show that every derivative of total order (time + space) up to n —2 of T vanishes on 
the boundary. The asserted regularity z G 5^ implies that did^z has zero trace for any j < n — 1 
and any k. Since any (n — 2)-order (space and time together) derivative of F involves at most 


(n — l)-st order terms in z, then (19) readily follows. Thus 

d{F G Cr^iA"^). 

Because ^4®/^ is by definition a bounded operator on ^(^®/^), then for | we have 

A^/^d{F = d{A^I‘^F G Ct^(A^^^) . 


confirming (14). 


□ 


Relying on (a special case of) Proposition 3.1 we formulate the notion of solution to 0-i)- 


Definition 3.4 (Weak solution to ([^-([^). We say u is a weak solution to Q -(§ on [0, T] if it is 
a weak solution to [[T, no,ni| with F{t,x) = f{u{t,x))u{t,x) (which is in L^LF‘{VL) by Proposition 


3.1 for n = s = j = 0). 


Definition 3.5 (Energy). For a function z = z{t,x) define the quadratic energy functional of order 
n by 

E^J^\t)-.= ]^\dfz{t)\\ + ]^\dfz{t)\l with (20) 
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3.2 Main theorems 


With the above definitions in mind, the new results of interest are: 

Theorem 3.1 (Global well-posedness of weak and regular solutions). Let f be as in @ 

with exponent 2m, m G N. Suppose for some integer n > 0 the initial condition Q satisfies 

(uo,Ui) G := ^ ^ 


Then there exists a unique weak solution u, regular of order n (Definition 3.3), of (§-(§ on [o,r] 


for any T > 0. Moreover u E SJf, that is, u E 2 ^) for j = 0,1,..., n. 


Excimpl© 3«2. If a solution to 0-(© with f{u) = u^ii has initial data given by the eigenfunctions 
of A (every ^(A’’)), then u G C'°°([0, oo); for every p > 1. 


Theorem 3.2 (Lack of uniform stability). The dynamical system generated by 0-@ on 
the state space corresponding to weak solutions := Hq{LI) x L^(fl) is non-accretive, but is 
not uniformly stable. Specifically, the energy functional t i—)• Eu{t) is continuous non-increasing; 
however, for any constants 0 < r <r and any time T > 0 there exists an initial datum {uo,ui) G 
Br{0) such that the corresponding solution trajeetory on [0,r] does not interseet 5^(0). 

However, for any s G (0,1) the lower-order norm ||n(t)||y(/s, 2 (Q) decays to zero os t —>■ oo with 
the decay rate uniform with respect to bounded sets of initial data. In addition, Eu(t) is strictly 
monotonically decreasing on any interval where Eu{t) is positive. 

Remark 3.1 (Open problem). The question of strong stability of the dynamical system generated 
by 0-0 on the finite energy spaee remains open, as far as we are aware. That is, for any 
given solution can we prove that lim^^oo Eu{t) = 0? 


4 Well-posedness 

4.1 Linear problem 

For convenience let us summarize a few classical results that can be easily verified, for example, 
using separation of variables. 

Lemma 4.1. Consider the problem [[F = 0, uo,mi| with {uo,ui) G Then there exist a group 
1 1 —)• S{t) of linear operators on such that {u{t),u{t)) := S{t){uo,ui) determines a weak solution 
to this problem on every [0, T], T > 0. The group is explicitly given by 


S{t){uo,ui) := '^2 


k=l 


cos('\/A^t) ^ sin(V^t) 

k sin (\/Afci) cos{^/Xft) 



Ek [uo] 


Ek 


_J^k[ui]_ 


_Ek_ 


( 21 ) 


where {{Xk,Ek)} are the eigenvalue-eigenfunction pairs for the operator A defined in 0, and Ek 
is the k-th Fourier coefficient with respect to the Hilbert basis {E^} for ififil): 

Ek{x) = V2sui{k'Kx) and \k = . 

Moreover, S{t) is a unitary operator on the space with respect to the norm given by as 

follows by direct calculation using (21). 
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The next result is likewise well-known: 

Proposition 4.1 (Inhomogeneous linear problem). Consider the problem [[P, rro, with 


{wQ,wi) £ . 


( 22 ) 


If F £ L‘^L‘^{Q), this problem possesses a unique weak solution w £ 


cO 

kD'j^ . 


Then the continuous 


mapping 1 1 —)• Eyj{t) for energy functional (20) satisfies 


Ew{t) + / {F{s),w{s))Qds = Ew{0) for all t£[0,r]. (23) 

Jo 

In particular, from the Gronwall estimate it readily follows 

\\{w{t),w{t))\\l^o < c (||(wo,u:i)||^o + ll-?^lli^L2(n)) e* for all t £ [0,r] . 

Moreover, if {wq,wi) £ and F £ then w £ and Wxx £ CtL'^{Ci) (e.g. see 

Thm. 2.1, p. 229^^ 


□ 


4.2 Variational formulation 

Proposition 4.2 (Variational form). Suppose u is a weak solution of ([^-([^ on [0,T]. Then for 
any test-function v £ Cll?‘{Q) n CtIlQ{Cl) with t £ [0,T], it satisfies the variational identity 


(h,'y)o 


/ («>^)o+ / / {f{u)u,v)^=^ 

Jo Jo Jo 


(24) 


Proof. Let c £ C^([0, T]), then for any 4> £ LIg(f2) we have from (12) 


(n(t),c(t)0)o 


{u{t),c'{t)(j))odt+ {u{t),c{t)(f)^dt+ {f{u{t))u{t),c{t)(j))dt = 0. (25) 


Let {En} be the orthonormal basis for lP‘{fiT) consisting of the eigenfunctions of A. Given v £ 
Cll?‘{il) n C'tifg(II) we can represent it as 


OO 

v{s,x) = ^Cfc(s)efc(x) 
k=l 


for Cfc £ (^^([OjT]), A: £ N. This series that converges to v in CtHQifiT) and its time-derivative 
YlT=i o'k(.^)ok{x) converges to v' in CtL‘^{0.). 

By applying identity (25) to finite sums ^ and passing to the limit M —)• oo we recover 

(M. □ 


^There’s a minor misprint in |24l Thm. 2.1, p. 229]: the first assumption is meant to read ug € Hglfi) (instead 
of “HglQ)"). In the second half of that theorem, which is the one we cite, it is strengthened to ug G Hg{Q,) n H^{Q) 
for strong solutions. 
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4.3 Relation between regularity and higher-order energy 

The existence of finite energy solutions to 0-® is known. Well-posedness for regular solutions, 
however, requires more work and relies on the connection between the smoothness in space and 
smoothness in time as summarized by the diagram below: 


,..,2 1, 


local solutions in 

t 

bound on the norm 


Et\t) := l\dMt)\'i + ■ 

—)• E^\t) is well defined 


I 


Energy identity and a bound on (t) 


The purpose of this subsection is to furnish this connection which can be loosely outlined as 
follows: a weak solution u of 0-0 is regular of order n on [0,r] if and only if the n-th order 
energy is bounded on [0,T]. We split this claim into two propositions. 

In order not to keep track of how the structure of f{u)u changes after differentiation we introduce 
the following, somewhat abstract property: 

Definition 4.1 (Regularity dependence). We say two functions {z,F) have order n dependence 
if for every k = 1,2,.. .n the regularity (if it holds) 

z G for j <k 

would imply that 

F G for j <k 

where j is a non-negative integer. It is helpful to note: 

• if{z, F) have order n dependence, they trivially have order i dependence for any i = 0,1,2,... ,n. 

• if (z, F) have order n dependence, then for s G N, s < re the functions A^l‘^z and A^^'^F have 
order (re — s) dependence. 


Again, it is helpful to consider an example. 


Example 4.1. Let F = f{z)z. Then it is not hard to check that {z,F) have order re dependence. 
For instance take re = 4 and assume z G ). As was shown via © and ( |18[ ) in Example 

T71 

d‘(F G &{A) 


continuously in time. In particular, F G C^^{A 2 ). 


Proposition 4.3. Suppose that u is a weak solution on [0,r] to [[F, reo,rei| with F G CtL‘^{LI) 
and \\{u,u)\\cj,je’o < Ri. Assume functions {u,F) have order re dependence for some re G NU {0}. 
Suppose, in addition, 

{dfu, dfu) G (26) 

with \\{dfu, dOu)\\c.j,^o < R 2 , then 

uGSlf, (27) 
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in other words, u G C^^{A ^2 ^ ) for j < n + 1. Moreover, 


iu,u)\\sji < K{n,Ri,R2), 


( 28 ) 


for a constant K{n, Ri, R2) dependent only on n, Ri, R2 and being continuous monotone increasing 
with respect to the parameters Ri, R2. 

Proof. Throughout the argument below, the norms in the considered spaces can be (inductively) 
estimated in terms of Ri and R2 , thus ultimately verifying ( 28). We will focus in detail on proving 
the claimed regularity. 

Cases n = 0, n = 1. By assumption we always have {u, u) G and F G CtL‘^{Q) which 

takes care of the case n = 0. If we in addition assume (ft, it) G then the equation il—Uxx = F 

implies that Uxx £ CtL‘^{^). Since tt = 0 on the boundary then u G Ct^{A). Conclude: 

(u, it) G . 

For n > 2 proceed by induction: suppose the result of this Proposition holds for 0,1,... ,n — 1. 
Assume (26) holds. Then let us show (p7|). 


Case n > 2. Because (rt, it) G then condition (26) implies (tt, it) G As a special 

case (using n — 1 instead of maximal n), we have 

Moreover, functions {u,F) have a fortiori (n — l)-st order dependence, so the induction hypothesis 
gives 

u G . 


Next, by assumption (26), we also have 

dfii = G CtL^{Pl) . 


(29) 


To show that u G Stf, it remains to verify that for j = 0,..., n we have dju ^ Ct^{A ^2 ^To 
this end introduce 

^ and F = A^/'^F. 

Since (tt, F) h ave n-th order dependence, then 2 : and F have (n — l)-st order dependence (see 
Definition |4.l[). As we already know, u G and via the (n — l)-st order dependence of (tt, F), 


we 


have F G C^^{A 2 ^ ) for j < n — 1. Because n > 2, then A^/'^F G likewise n > 2 

also tells us that u is at least in Ct&{A), so AAI'^Uxx = A^^‘^u = dxxA^^'^u G Ct&{A~^^‘^). 

Consequently, applying AA!’^ to the equation for u, we conclude that z = AA^'^u is a weak 
solution to [[T, 2 ( 0 ), i;(0)| with F G CtLA{^) and ( 2 ( 0 ), i:(0)) G By (26) we also have 

The induction hypothesis now states that 


z G a: 


n—1 


for j = 0, l,...,n. 


It implies that 


^Cj.S>{A ^2 for j = 0,I,...,n, 


u 


as desired. From here, along with (29), it follows that u £ Sf, which completes the proof of the 
implication “case n — 1 ^ case n” for the induction argument. □ 
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The next result complements the previous proposition, demonstrating that the same regularity 
of solutions can be inferred from a priori differentiability in space, rather than in time. 

Proposition 4.4. Suppose u is a weak solution on [0, T] of [[P, uq, ni| with F G CtL‘^{^)- Assume 
{u,F) have n-th order dependence (Definition f.l) for some n G N U {0}. If u is regular of order 
n, i.e., \\{u^u)\\cj.M’‘^ < R < oo, then 

u&SJf 


and 


(30) 


||('w,'u)||s?; < K{n,R) 

with K dependent only on n and R and continuous monotone increasing with respect to parameter 

R. 


Proof. In the course of the proof, the bound (30) can be traced inductively to depend only on 
\\{u,u)\\cj,j)f". To keep the exposition concise the argument will focus on the regularity. 

If re = 0, the claimed regularity simply matches that of weak solutions. For re = 1 we are 
given u G CtS’{A) and u G Ct^{A^/'^). Solving u = Uxx — F £ verifies that u G 

CT^{A)nC^L^{n) = S^. 

Proceed by induction. Fix re > 2, suppose the statement holds for all re < re — 1 and assume 


II(re, u)\\cj.jffri is finite. A fortiori it follows that (re, re) G Then re G by the induction 

hypothesis. 

Define z := A^/^re, then using re G it follows as in the proof of Proposition 


4.3 


that z 

is a weak solution to z(0), z(0)| with Afil’^F G CtL‘^{VL). Moreover from the assumption 

re G we also have that z G i.e., is regular of order re — 1. Thus by the induction 

hypothesis 


z G 5 ; 


n—l 


equivalently 


£ C^&{A 2 ) for j = 0,l,...,re. 


re 


(31) 


The only remaining step from here to proving re G 5^ is to show that £ CtL‘^{VL). To this 

end define 


w 


■.= d. 




u. 


Then by (31) 


(re;, w) £ ■ 


Since u £ SJf ^ and by the (re — l)-st order dependence (implied by re-th order dependence) of re 


and F, we deduce from (31) that df ^F £ CtL^{D). So 

-d'^-^F £ CtL^D) 


W = Wx 


confirming that 

w = £ CtL^{D) 

as desired. Thus u £ SIf completing the induction argument. 


□ 
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4.4 Local unique solutions 

It is known that 0-0 possesses unique solutions [l8]. Here we extend this result to regular 
solutions as well. First we formulate it for local in time solutions. 


Theorem 4.1. Suppose 11(^0)1^1)11.^" = R < oo for a non-negative integer n. Then there exists 
T = T{R) > 0 such that system 0-0 has a local unique solution that is regular of order n on 
interval [0, T]. 

Proof. Step 1: the spaces. Note that for if G 

IV^Io + Idxi’lo ■ 


< 

rsj 


(32) 


Moreover, for (ijjQ,ifi) G we have 

with the temporary notational convention (^“^(H) := L^(H). 

Step 2: contraction mapping. Let t S{t) be the semigroup for the linear wave equation 
with data go G Note that if z is regular of order n, then by Proposition 3.1 we would get 


f{z)z G . 

With this in mind introduce the operator At on 

0 

- s) f 

/o 

For t G [0, T] we have 


(33) 


{^T(f){t) ■■=S{t)yo+ [ S{t-s) 
Jo 


fiMs))Ms) 


ds for any f = (^q, fi) G CtA^^' 


AT(j){t) — AT^{t) 


./o 


:/ 


S{t — s) 




fiMs))4>i{s) - fiMs))4>i{s 
From the local Lipschitz property ([^ of / we obtain 

f{(fo)(pi - /((^o)0i = f{(fo)(fi - f{^o)4'i + f{^o)fi - fiMfi 
=cj)iM{4>o,(j)o){(t)o - ^o) + f{^o)i(l)i - 4>i) . 

Now we will rely on the fact that a priori </>,</> G gives: 

00, 00 e CTC"(n), 01, 01 G and 0i,0i G CrW^'^in). 

Estimate, 

/(0O)01 - /(0O)01 


ds. 


(34) 


^(A"/2) 


< 


< 


0lM(0O,0o)(0O - 0o) + /(0o)(01 - 0l) + 0”[0lM(0O,0o)(0O - 0o) + /(0o)(01 - 0l)] 
0lM(0o,0o)(0O - 0o) + /(0o)(01 - 0l) 


+ 

Cijk dl(j)i ■ ( 9 ^M( 0 o, 0 o) • 9^(00 - 0 o) 

+ 

Y ^b ' 9 x/( 0 o) • 5 ;J( 0 i - 0 i) 



0 


<Vi 

1 h/’O C'n(Qp 00 C'n(Q)| • 01 W"’2(n) • 00 ~ 

*^0 0" (fi) +"^2 1 00 c"(n)} ■ “ 


( 35 ) 
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where 'Pi{-} are a polynomials in the indicated variables with positive coefficients. Now, from (32) 




< Wfh- 




So continuing (35) we get 


fi4>o)(pi - /(<^o)(^i 


<vA f 11*^ “ 


I 


Because the group S{t) is unitary on ^(A"") then (34) yields 
h.j'cj) — A.j'cj) 


< 




TP 4 \\(t> - ■ 


(36) 


If (p and (j) come from a bounded set, then choosing T small enough implies that is a contraction. 
Step 3: invariance. Finally, we also need A.t to map the the ball Bji{0) in into itself. 


(37) 

□ 


According to (36) for that it suffices that R satisfy 

TVi{R, R} 2 R < R . 

For which we can impose 

T< l/(21P4{i?,i?}). 

Now the contraction mapping principle implies the claimed local unique solvability. 

4.5 Energy identities and global existence 

The global existence result is based on a priori bounds on the energy. 

Proposition 4.5 (Extension of local solutions). Let f be as in (Q with exponent 2m, m G N. 
Assume u is a weak solution that is regular of order n, defined on some right-maximal interval 
Imax = [0,Tmax)- Suppose {u,u) G and that there exists a continuous function ip on M+ 

such that 

Ejfi\t)<ip{t) on 

Imax ■ 

ThcTl TyrifL'K — 00. 


Proof. By Proposition 3.1, u and the term F = f{u)u have n-th order dependence. Recall that 
energy controls the norm of the solution (in fact, it controls the norm, but 

the continuity with values in is implied by the definition of regular solution). Then there is a 
constant K dependent on sup{|V’(t)| : 0 < t < Tmax} such that 

||(u,h)||crj^" < K 


for any T < Tmax- Hence by Theorem 4.1 any initial data of the form vo = u{t) and vi = iiif) 


for t G [0,Tinax) can be extended to a solution that exists for another AT = AT{K) time units, 
independently of t G [OjTmax)- Hence T^ax cannot be finite. □ 
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Lemma 4.2 (Energy identity for regular solutions). Let f be as in Q with exponent 2m, m G N. 
Suppose u is a weak solution to on [0, T]. If u is regular of order n, then for k < n we have 

{u,u) G f{u)u G and 


£<'■'(()+ / (8‘[/Wi|,8f!i),=EW(0) 


(38) 


Proof. Proposition 3.1 confirms that F and F = f{u)u have n-th order regularity dependence. 
Invoke Proposition 4.3 to conclude that u ^ Stf. It certainly implies that (u, u) G for k < n. 


And by the re-th order regularity dependence of {u, F) we also have F G C^L‘^{Q). 
Define 


w 


= dU 


then ui is a weak solution to \d^F, r(;(0), u;(0)| with d^F G CtL‘^{Q) and (t(;(0), t(;(0)) G Now 
call upon the energy identity (l2^ for w to obtain (38) for u. □ 


4.5.1 Finishing the proof: global existence and regularity 

Let tt be a weak solution on interval [0, T]. Because F ■ ii ■= f{u)vf is non-negative by assumption 
@ on /, then in the case n = 0 by Lemma [4.2| we have 


\\\{u,u)\\l; ^0 < sup Eu{t)<Eu{^). 
te[o,r] 


(39) 


Since the bound is independent of T G [OjTmax), then by Proposition 4.5 u extends globally to 
t > 0. 

Arguing by induction, suppose that (u, ft) is a global solution that is regular of order re, and at 
the same time is local regular of order re + 1 on a maximal interval [0, Tmax)- We want to show that 
it is also global of order re + 1. 

Example 4.2. It helps to preview the result on the example of f{s) = and by extending from 
weak to regular of order 1 solutions. Suppose u is a global weak solution that is, in addition, regular 
of order 1 on maximal interval [0, Tmax)- For such a solution by Proposition \4.i^ we have the “1-st 
order” energy identity 

Ei^\T)+ r dt[f{u)u]u = Ei^\o). 

Jo 

rT 

eL^Kt)+ / {2uu‘^ + u^u)u = Ei^\o). 

Jo 

We can bound CtL°°{I}) norms ofu and ii in terms of the finite energy Eu(t). The latter is bounded 
by constant Eu{0) (but, in fact, any continuous on M+ upper bound on Eu{t) would do, which is 
how the general argument works). In turn, the term |re|^ is a part of the first-order energy. So we 
have 

E(i)(r) < eW( 0) + C(Ef(0)) r Ei^\t)dt 

Jo 

From here, Gronwall’s es tima te gives us an asymptotically growing continuous upper bound on 
E^\t). Then Proposition f.S ensures that u is regular of order 1 globally, that is Tmax = oo. 
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Now onto the actual proof of global existence. By the n-th order global regularity there is a 
continuous monotone increasing function ipn ■ IR+ —>■ K+ 

\\{u,u)\\sr^ < for all t>0. 


By (n + l)’st order regularity and via Lemma 4.2 for any T < [0, Tmax) we also have 


4 ^+ 1 ) (T) + / {dr^F, dr^u)o = 

Jo 


Use estimate (15) of Proposition |3.1| to claim that there is a polynomial V{s) such that for each t 

lar'Fwi. < rmmi +isr‘«(*)io) ■ 

Then |ar+‘F(t)U8r+‘i.|. < P(V-„(*))(1 + < 2P(V-„(*))(1 + |8r‘«i;). In 

particular, 

4-+i)(T) <4"+i)(0) + 4 rV{^Pn{tm + Et^"\t))dt. 

Jo 

Gronwall’s inequality now provides 

< i^n+i{T) := + j\v{Mt))d^ exp (^j\v{Mt))dt^ for all T < r„,ax 

By the assumption on the newly defined mapping "^n+i is continuous on M_|_. Then Proposition 


4.5 ensures that the solution u is global, regular of order n + 1, which completes the proof by 


induction. 


This step also completes the proof of Theorem 3.1 The additional regularity u G follows 


if we first use Proposition 3.1, which establishes n-th order dependence of u and f{u)u, and then 


invoke Proposition 4.3 which translates n-th order regularity into S^. 


5 Lack of uniform stability 

This section is devoted to the proof of Theorem |3.2[ The heart of the argument is that the energy 
decay slows down with the frequency of solutions, independently of their total finite energy. 

Let k and consider initial condition of the form 


no^^(0,x) = sin(A;7rx), dtu\^'{0,x) = 0. 
kiT 


.{k)i 


(40) 


Note that u'^'^ G &{A^) for any r G M. By direct calculation we have 


= 1 and \\uf\\ = ^. (41) 

Hence for such initial data, any constant that may be estimated in terms of ||(ug^\u^^^)||jfo is 
bounded uniformly in k. 
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5.1 “Primitive” problem and decay of lower-order norms 

,(fc) 


We continue working with the initial data (uq ,0) as in (40). Pick the anti-derivative through the 
origin of / (as in Q) 


Rs) := 


„2m+l 


2m + 1 


Due to IRRiR' (hence being bounded independently of k, we have for some constant 

C independent of A: in (40) the estimate 

\f{Ro\<C\Ro^\o- 

Let function solve the linear elliptic boundary value problem 


f = f{RR) 

I $(fc)(0) = 4>W(1) = 0. 

It satisfies the standard elliptic estimate: 


(42) 


(43) 


for C' and C independent of k. We do have G L^(D) because RR G H^y). In fact, u, 

is an eigenfunction, and it is easy to show (as in Proposition 3.1) that 


(fc) 

0 


7(4"^) G^(A)nC'~(D). 


We then have 


yk) ^ ^ W^’2(D) . 

With this function at hand, consider the nonlinear “primitive” problem 

— (jR} + 7(7*'*'^) =0 for X € Q, t > 0 
with homogeneous Dirichlet boundary data 

= = 0 for t>0 

and initial condition 

y\o) = and = uo ■ 


(44) 

(45) 

(46) 


Note the following properties of the system (44)-(46): 


i.) Because s i—?• /(s) is continuous, monotone increasing vanishing at 0, then (44) is a semilinear 
wave problem with monotone damping. The data (<I>^^74^^) come from the domain of the 
associated nonlinear generator; in this ID setting it coincides with the domain of the 
corresponding linear group. This well-posedness result is based on monotone operator theory 
[251 Thm. 3.7, p. 306]. 

In particular, we have for any T > 0 

^(k) ^ n and G LfL^Q). 


19 









Moreover, given the ID embeddings it follows that 

G C^L'^iQ) n . 

So satisfies (44) and has ((/>(0), (^(0)) G ^{A^) x By Proposition 

G CtL^q) . 


4.1 


we have 


ii.) The velocity v^^\t,x) := satisfies G and is a weak solution to 


vW - vi’^J + ^ = 0 


In addition, 


tj(*^)(0) = (^•^^^(0) = Uq^'^ G L^{Vt) (actually, in any S>{A')) 


(47) 


and because are continuous with values in L^(D), then 


•,(D 


(0) = <)iW(o)® 4"i(o)-/(<^('^)) = a2ck(")-/(4'=)) = o 


according to the way 4>o was constructed in (42). Thus the primitive problem is the velocity 
potential for our original problem. Namely the time derivative is the solution to our 
original degenerately damped system with initial data (uq^\o). 


hi.) Solutions to (44)-(46) decay uniformly to zero at the rate that can be estimated explicitly 
in terms of the exponent 2m of /. This is a fundamental example of a dissipative problem 
with full interior damping. The decay rate can be assessed using, for example, the ODE 
characterization of [22] (see jS] Coro. 1, p. 1770] for more details; in particular function h(s) 
there must satisfy /i(s^™+^) > s^). We get, as t —)■ oo that 

1/m 

E^(k) it) 


It follows that the (lower-order) norm |M^^^(t)|p = decays to zero as t —>■ oo with the 

decay rate uniform with respect to the finite energy E^(k) (0). By interpolation between L^(D) 
and we conclude uniform a decay (but at slower rates, where 1/m is modified by the 

interpolation exponent) for every norm with s G (0,1). 

Remark 5.1. Note that uniform stability of the original system would have required us to 
include the case s = 1, which as we are about to show cannot happen. 


Summarizing: is precisely the solution to ^f 0^. Any fractional 

IT®’^(D) Sobolev norm of for s G (0,1) decays uniformly with respect to E.^(k){0). In addition, 
we have the estimate 

(| 43 |) 

|^(D(i)|2 = |</(fc)(t)|2 < 2E^,k,{t) < 2E^(.)(0) = V C2\ui/^\l (48) 

for C 2 independent of hence independent of k in ( |40| ). In particular, the lower order norm 

|rt(*^)(t)|Q of the solution to our original problem ^f {u)u^^\ is controlled by its initial 

(fc)l 

lower-order norm \uq independently of k. 

Now we are going to use the smallness of the L^(D) norm of A^\t) to show that its ^(A^/^) 
norm cannot decay too fast. 
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5.2 Comparison with a conservative problem 


Let be the unique solution to linear homogeneous problem [[0 ,Mq^\o|. From the energy 
identity (23) and by the choice of in (|40j) we get 




Define 


Zik) = uW _ yj(k) 


then solves The energy identity ([2^ for 2 ; gives 




=0 


Note that 

• ( 0 )-E^{fc) ( 0 )) < 2 C (since -E„{fc)( 0 ) = 1 ) independently of k in 

• Given / as in Q we have Via 1-dimensional embeddings and 

interpolation for es < 1/2 


\u 


(fc)| 


L°°(n) 


< 


,(D| 


wi-=.2(o) 


< 


u 


(k) 


u 


(k) 1 1—£ 


SO 


< C(E„(«(0),e) k'*’!' 


Because Eu{0) = 1 independently of /c in (40), then by (48) we get for any t > 0 

ll/(«“>(i))llL»(n) < < C3|«i 


,{k)\£ 
^0 lo- 


for C 3 independent of k in ( |40[ ). At this point we finally expand the definition of uq to get 
(see dil])) 

\\f{u'^^Hi))\\L-°{n)<C 3 ,£^ for any e < 1/2 and f> 0 . 


Plugging these observations into the estimate for E^(k) (t) we obtain: 

Lemma 5.1. For u'^\x) '■= 2 (/c 7 r)“^ sin(fc 7 rx) let denote the weak solution on [0,r] to the 
problem and be the solution to linear homogeneous problem [[ 0 ,rig^\o|. 

Then the differenee — rcV) solves ,0, 0| and for any e <1/2 it satisfies 


E^(k) (t) < CsT 


1 


for all t G [0, T] , 


with Cs independent ofk. 


□ 
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5.3 Ruling out uniform stability 


At this point for brevity let us suppress the superscript “(fc)’’ 
parameter in the dehnition of initial data. Let u,w,z be as in Lemma 5.1 
that Ez{t) < 5 for some t. Then via E^{t) = 1 = Eu{0) we have 


stemming from the choice of the 
Suppose for a moment 


Eu{t) = + Zx{t)\l 

= Ew{t) + Ez{t) + Zt{t))^ + {Wx{t), Zx{t))^^ (49) 

> E^{t) + £' 2 (t) - 4y/E^{t)Ez{t) > 1 - 4VS . 


To apply this estimate, pick any T > 0 and hnd 5 > 0 such that l—4y/6 >1/2 (e.g., if J < 1/64). Fix 
s < 1 / 2 , then there is k = k{S) large enough so that for initial condition {uq = 2 ( 7 rA:)“^ sin(/c 7 rx), 0) 
yields a solution whose energy satishes 


E,{t)<Ce^<6 for fG[0 ,r], 


for Ce as in Lemma 5.1 Consequently, by (49), 

Eu{t) > - for all t G [0, T] . 

However, the initial condition {uq,0) had energy 1 (again, independently of k). Hence the family 
of initial conditions 

(50) 


I := ^ sin(fe 7 rx), '■ k G N 


with the associated solutions resides in a bounded ball (of radius ■s/2E^(k) ( 0 ) = v/ 2 ) in Jif°, 
yet the corresponding solutions do not decay to zero uniformly in the topology of 

Thus the associated dynamical system on is not uniformly stable. This argument demon¬ 
strates Theorem 3.2 for f = 1 and r = ^. The general case follows merely by attaching a factor of 
,(fc) 


r to Uq and choosing a potentially smaller 6 in the last step of the argument. 


5.4 Monotonic strong decay of the energy 

For a weak solution u of ^f{u)u,uo,ui^ with (rio,ui) G the functional 

1 1 —>■ Eu{t) 


is non-increasing. We cannot presently claim whether solutions decay to 0 strongly, however, it is 
possible to show that the energy Eu{t) is strictly monotonically decreasing for non-trivial solutions. 

Consider a weak solution u on [0, T] and suppose Eu{t) = const on (a, 6 ) C [0, T], then from 
the energy identity (23) follows that 



it))^dt = 0 


Thus u’^ii = 0 a.e. in (a, b) x H. Since it is equivalent to = 0, then we conclude that 

ri™’+^(ti) = ri™'+^(t 2 ) in L?{Q) for a.e. ti,t 2 '. 


= const in L^(H) a.e. t G {a,b). 
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Moreover, since u G CtHq{Q) “—)• (^([OjT] x fi), then we have u{t,x) = ±no(x) for every {t,x) G 
(a, b) X O. By the continuity of the solution this is only possible if u = uq for (t, x) G (a, b) x fl. 
Then u = 0 and we arrive at an equilibrium solution which has to be trivial. This observation 


completes the proof of Theorem 3.2 


6 Numerical results 

The strategy for the proof of instability was largely prompted by numerical observations described 
below. 


6.1 Outline of the numerical approach 

The numerical implementation presented here treats the case m = 1 of Q: 

o 

u — Uxx + u u = Q 


with 

and given initial data 


u{t, 0) = u{t, 1) = 0 

u{0,x) = uo{x) and u{0,x) = ui{x). 


Solution was discretized in space via a Ritz-Galerkin method. The dynamic problem could be 
analyzed explicitly using a discretization in time and a Runge-Kutta scheme, though, rigorous jus¬ 
tification of convergence becomes more delicate. Another approach is to approximate the successive 
approximations of Theorem 4.1 which, when exact, are guaranteed to converge, at least over small 
time intervals. The iterates correspond to linear inhomogeneous PDE problems that are resolved 
using a hybrid scheme: 


(i) for relatively short times find solutions using an approximation of semi-discrete Ritz-Galerkin 
method by discretizing time-integrals in the variation of parameter formula. If the error in 
numerical integration is small, then this approach enjoys an explicit convergence estimate 
(for smooth solutions and over finite time intervals) essentially proportional to the space 
discretization parameter h. 


(ii) for larger times, collect the last /c-points of the semi-discrete approximation and resolve the 
rest of the iteration using a multi-step method (Adams-Bashforth). 


Thus, we begin with some initial guess (up],ujo]) and proceed to solve inhomogeneous linear 
problems 

^[k] ~ dxxU]j~'^ = ( 51 ) 

on the space = Hq{Q) x with the forcing term from the preceding iteration 


F[k-i] - 


= —Ur: 




The constants in the estimates (36) and (37) could potentially be determined explicitly in this 
case (by following the proof with specific a and m in the definition (Q of /) which would yield an 
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explicit bound on the Lipschitz constant “ 7 ” of the contraction mapping in terms of T. In turn, 
given this constant 7 > 0 if 

\\u[k] - U[k-i]\\cTJ^° < ^ 


then the absolute error between k-th. iteration and the true solution is no more than £ 7/(1 — 7)- 
If we use a Ritz-Galerkin scheme with element size h to find an approximate solution U[k,h] to 
linear inhomogeneous problem (51), then for h small, e.g. see |26l Thm. 13.1, p. 202], and for 


simplicity taking the initial conditions to be the more accurate Ritz projections of the initial data, 
we get 


\u[k] - uik,h]\\cTJlf0 <Ch \u[k]{s)\ 2 ds. 

Jo 


This estimate of course requires sufficiently regular solutions. As Theorem 3.1 and Example 3.2 


show, in order to have the L^IT^’^(r 2 ) regularity on it suffices to have initial data in the space 


(u(0),u(0)) e ^ ^(A^) X . 


For example, the demonstrated numerical results below use displacement and velocity proportional 
to the eigenfunctions of A, which are smooth. 


6.2 Specifics of the implementation 

Consider an equipartitioned mesh of subintervals of length h and the standard piecewise linear 
nodal basis with i = := (h~^ — 1). Let A/j denote the restriction of the linear 

evolution generator A to the subspace of Jf" spanned by By TZ\ denote 

the elliptic Ritz projection on the subspace of Hq{Q) and let TZ^ stand for the corresponding L^(n) 
projection. 

Given a candidate approximation 


y[k—l,h] 


^[k—l,h] 


we compute the forcing 


f[k-l,h]{t,x) = uf^_^^{t,x)U[k-l,h]{t,x) 

The coefficient vector i[k^h] of the projection of Tl^fik^h] is obtained in terms of the coefficients 
Cpqrs = [ 4>^ix)^’^{x)(j)^{x)(j)'^,{x)dx (which for this choice of basis functions (f)'} form a very sparse 

Ju 

tensor with only 3 distinct values). The initial guess used to calculate /[o,/i] is the constant solution 


f[ 0 M^'^^^) = ^'^hUo{x)f{nlui{x)) for all t > 0 . 

Let yo denote the projection of the initial data {JZ\uo,TZ^j^ui). We obtain a semi-discrete approxi¬ 
mation of the original system for unknown coefficient vector y[k,h]- 

y[k,h] = -^hy[k,h] + f[k-i,h], y(0) = yo • 

For relatively short times we can invoke 

yk,h{i) = e~^^^yo + It with h ■= f e~^^‘'*'~'''>iik-i,h]{s)ds 

Jo 
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which is in turn discretized over time scale T = ■ ■ ■, tj\f) with ti+i —ti = 6. According to 

It = 

at each step only the integral over [t, t] needs to be computed. For this purpose only several 
values of the matrix exponentials are needed in order to apply the Newton-Cotes rule on 

sub-interval specifically: 

gi-<5Ah £qj. j = 1^ 2,..., m — 1 

where m is number of points used for Newton-Cotes formula (e.g., Boole’s or Simpson’s 3/8th). 
These m — 1 matrices need to be computed just once and only depend on the time-step, but not 
the total number of these steps. 

In turn, the vectors have to be found for each But since yo is fixed, these can be 

more accurately determined using scaling and truncated Taylor series approximation m- 
For simulations over larger time-scales we can use the last few values: 

{tN-p, yk,h{tN-p)), ■ ■ ■ , {tN, yk,h{tN)) 

to initialize a linear p-step method, e.g., 5-step Adams-Bashforth to efficiently obtain solution on 
the interval [tAr,ifinal]- 


6.3 Pointwise Runge-Kutta solutions for particular initial data 

As before, let be the eigenfunction \/2 sin(7rA:x) for A with eigenvalue \k = Then for 

initial data 

uq{x) = coEk{x) and ui = ciEk{x) (52) 

for constants cq, ci, the solution of 0-@ can be reduced to a dissipative ODE using the ansatz 

u{t,x) = (i){t)Ek{x). (53) 

Plugging it into 0 equation yields 

(t>"llk + ^k4>IIk + Ek4>^ 4>' = 0 


This identity would be implied if for each x G D function (f> solves the 2nd-order nonlinear ODE 

(/>" Afc0 -h Ek{x)(j)‘^(j)' = 0 . 


i?!)(0) = Cl and 0^(0) = C 2 . 
It corresponds to a first-order nonlinear system: 


T' = Ffc(T,x), 


f 

a 

\ . 

/3 

[ 

P_ 

,x 1 : = 

-Xka - Ek{x)a^/3_ 


T(0) 


Co 

Cl 


(54) 


Function F^. is smooth with respect to the components of T and to variable x, which now acts as a 
parameter. This ODE system has global differentiable solutions, moreover since Ek{x) is smooth, 
in fact, analytic in x then local solutions are differentiable with respect to x |28l Thm. 3.1, p. 95]. 
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Figure 1: The displacement value u{t,x = 0.5) of a numerical solution to problem Q-Q with 
f{s) = Initial displacement uo{x) = 27r“^ sin(7rx), initial velocity ui{x) = 0. Time interval: 
[0,r = 10]. Obtained by approximating the Ritz-Galerkin semi-discrete solution using numerical 
integration in the variation of parameter formula. Element size: h = 10“^; time-step: 6 = 2- 10“^. 
Displayed next to the exact analytic solution of the corresponding initial boundary value problem 
with linear damping: u — Uxx + (2/7r)^tt = 0. 


Because the initial data (52) is smooth then by Theorem 3.1 the unique solution u is, among 
other things, in CyC^(D). Consequently u must coincide with the solution to the ansatz (53). 


In turn, (54) is a dissipative 2x2 system of ODEs and can be approximated by a Runge-Kutta 


scheme. To get some quantitative estimate on the absolute error of solutions found in Section |6.2| 


at least for initial data of the form (52), one can consider a piecewise linear interpolation of 


and then calculate the energy-norm difference from the finite-element solution. 


6.4 Energy plots 

The accompanying figures and data demonstrate some of the numerical results. The initial data is 
considered of the form 

2 

uo{x) = — sinlfeTTx), ui = 0 
nk 

which permits to compare the finite element solutions to the pointwise Runge-Kutta solutions 
described in Section 16.31 

Eigure shows the point-value of displacement u{t,0.5) next to the displacement value at the 
same x = 0.5 for the corresponding initial boundary value problem with linear damping. 

Eigure presents numerical estimates of the energy for solutions obtained by Ritz-Galerkin 
finite element scheme and successive approximations. The graphs indicate that the energy decay 
deteriorates as the frequency of the initial data goes up while the initial finite-energy remains fixed 
(E^(fc)(0) = 1 independently of k), thus illustrating the lack of uniform which was rigorously con¬ 


firmed by Theorem 3.2 The initial data are of the form (52) with zero initial velocity. The indicated 


errors are obtained by comparing each finite-element solution to a piecewise-linear interpolant of 


the corresponding piecewise RK solution (53). 


Eigurej^uses multi-step extensions of the same solutions shown in Eigurej^to a larger time-scale 
using (5-step) Adams-Bashforth method. It also includes the decay of the L^(D)-norm \u^^\t)\g 
for these solutions. 
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Decay of the energy E(t) on [0,10] 



Figure 2: Plots of energy E^{k) {t) versus time t for system 0-@ with f{s) = Time interval: 
[0,r = 10]. Obtained by approximating the Ritz-Galerkin semi-discrete solution using numerical 
integration in the variation of parameter formula. Element size: h = 10“^; time-step: 6 = 2- 10“^. 
Initial displacements: = 2(A;7r)“^ sin(A:7ra:) for frequencies k = 1,2,4, 8. Initial velocities are 

zero. Approximate error was obtained by computing maximal, over [0,T], energy ( 0 ^*^-norm) 
difference from a piecewise-linear interpolation of the pointwise RK solution given by ansatz (53). 
The initial energy of each solution is 1. Corresponding errors are: ei = 3.84E-02, 62 = 5.20E-03, 
64 = 6.97E-03, 68 = 4.15E-01. 




Eigure 3: Plots of energy E^(k) (t) and norm vs time t for a numerical solution of 

problem ([^-([^ with /(s) = s^. This is an extension of the solutions from Eigure (originally 
defined for t G [0,10]) to the time-interval [0, 50] by 5-step Adams-Bashforth method. Due to rapid 
oscillations, the L^(D) norm values for the solution corresponding to A: = 8 appear to fill out a solid 
region. 


27 

















































Acknowledgments 

The research of George Avalos and Daniel Toundykov was partially supported by the National 
Science Foundation grant DMS-1211232. The numerical analysis presented in this work and partial 
theoretical results were obtained during a Research Experience for Undergraduates (REU) program 
on applied partial differential equations in summer 2013 at the University of Nebraska-Lincoln. This 
REU site was funded by the National Science Foundation under grant NSF DMS-1263132. The 
REU project would not have been possible without assistance of Tom Clark. 


A Finite-dimensional counterpart 

To complement the analysis of the infinite-dimensional model Q-Q it is also interesting to examine 
the related finite-dimensional version of a degenerately damped harmonic oscillator: 


X + kx + f{x)x = 0, {k > 0) 
x(0) = xq, i:(0) = xi 

for / = m G N. Rewrite it as a first order evolution problem 


y = G{y) for G{y) := 


2/2 

-f{yi)y2 - kyi 


, 2 /( 0 ) = 


XQ 

Xl 


Henceforth, let | • | denote an equivalent norm on 


lup := \k\vx\^ + l\v2\‘^ . 


(55) 

(56) 

(57) 

(58) 


Because / is smooth and non-negative, then classical ODE results guarantee that solutions are 
unique, exist globally and satisfy 


|2/(i)| < |2/(0)| for all t>0. 


Below we present stability results which contrast their infinite-dimensional analogues discussed 
earlier. In particular, the finite-dimensional system is uniformly stable, while for distributed- 
parameter version the strong stability is open while uniform stability has been proven false. 


Lemma A.l. The dynamical system corresponding to the ODE (55)-(56) is asymptotically (“strongly”) 
stable. 


Proof. This is a direct consequence of LaSalle’s invariance principle with the Lyapunov function 
given by the equivalent norm (58): V{y) = \y\‘^. We only need to check what kinds of trajectories 
reside in the invariant set of the system: 

E = {y€ 


: VV{y) • G{y) = 0} 


So for y G E 


kym - f{yi)yl - %i2/2 = 0 . 
f{yi)yl = 0 
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In particular, either yi = 0 or ?/2 = 0. 

Now suppose that a solution trajectory {{t,y = {x{t),x{t))) : t > 0} resides in E. If at some t 
we have x{t) ^ 0, then by the continuity in time, it is nonzero on some interval I. On that interval 


we must have i; = 0 by the property of E. But then on that interval, from equation (55) we get a 
contradiction since the solution has to be constant and therefore zero. 

Thus the only trajectory in E is the trivial one. By LaSalle’s invariance principle every solution 
is asymptotically stable. □ 


Theorem A.l. The dynamical system corresponding to the ODE (55)-(56) is uniformly stable. 


Proof. Proceed by contradiction. Assume for a bounded set B C there exists some e > 0, a 
bounded sequence of initial data {yon)n C B, and a sequence of corresponding times {Tn)n with 
Tn oo such that 

\yn{Tn)\ > S. 

Extract a convergent subsequence of initial data, reindexed again by n. Let zq & B denote this 


limit point and 1 1 —>■ z{t) be the corresponding solution. By Lemma A.l 


lim z(t) = 0 . 

t^OO 


In particular, there exists T such that for t > T we have \z{t)\ < e/2. Because the non-linearity 
(xi, X 2 ) e->■ f{xi)x 2 is locally Lipschitz on M^, and the system ( [ST] ) is non-accretive, then there exists 
6 = 6{T, B) > 0 such that for any rjo G B ii \ri — zo\ <5, the corresponding solutions satisfy 

\r]{t) - z{t)\ < ^ for tG[0,T]. 

Since T„ 00 , we can find T^r > T. Next, because j/on converge to zq, then for n > A we can hnd 
hnd yon so that |yon — zo\ <6 and, consequently, 

\ynit) - zit)\ < ^ for te[0,r]. 

Then \yn{T)\ < e/2 + | 2 :(r)| < e. Because the system is non-accretive, then |yn(0l < ^ 
t >T, and in particular it holds for T„ > Tj^ > T. So \yn{Tn)\ < e which contradicts the choice of 
yon- □ 
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